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Abstract 

The paper describes an extension of the Liga algorithm for structure solution from 
atomic pair distribution function (PDF), to handle periodic crystal structures with 
multiple elements in the unit cell. The procedure is performed in 2 separate steps - 
at first the Liga algorithm is used to find unit cell sites consistent with pair distances 
extracted from the experimental PDF. In the second step the assignment of atom 
species over cell sites is solved by minimizing the overlap of their empirical atomic 
radii. The procedure has been demonstrated on synchrotron x-ray PDF data from 
16 test samples. The structure solution was successful for 14 samples including cases 
with enlarged super cells. The algorithm success rate and the reasons for failed cases 
are discussed together with enhancements that should improve its convergence and 
usability. 
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1. Introduction 

Crystallographic methods of structure solution are the gold-standard for determining 
atomic arrangements in crystals including, in the absence of single crystals, structure 
solution from powder diffraction data (Pecharsky & Zavalij, 2005; David et al., 2002). 
Here we show that crystal structure solution is also possible from experimentally 
determined atomic pair distribution functions (PDF) using the Liga algorithm that 
was developed for nanostructure determination (Juhas et al., 2006; Juhas et al., 2008). 
The PDF is the Fourier transform of the properly normalized intensity data from an 
isotropically scattering sample such as a glass or a crystalline powder. It is increas- 
ingly used as a powerful way to study atomic structure in nanostructured materi- 
als (Billinge, 2008; Egami &; Billinge, 2003). Such nanostructures do not scatter with 
well defined Bragg peaks and are not amenable to crystallographic analysis (Billinge 
& Levin, 2007), but refinements of models to PDF data yield quantitatively reli- 
able structural information (Proffen & Billinge, 1999; Farrow et al, 2007; Tucker 
et al, 2007). Recently ab initio structure solution was demonstrated from PDF data 
of small elemental clusters (Juhas et al., 2006). Here we show that these methods can 
be extended to solve the structure of a range of crystalline materials. 

Whilst it is unlikely that this kind of structure solution will replace crystallographic 
methods for well ordered crystals, this work demonstrates both that structure solution 
from PDF data can be extended to compounds, and that robust structure solutions 
are possible from the experimentally determined PDFs of a wide range of materials. 
We also note that there may be an application for this approach when the space-group 
of the crystal is not known, as the Liga algorithm does not make use of such symmetry 
information. In fact, the space group can be determined afterwards analyzing the sym- 
metry of the solved electron density map (Palatinus &; van der Lee, 2008). However, 
this approach is promising for the case where the local structure deviates from the 
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average crystallographic structure, as has been observed in a number of complex crys- 
tals, for example the magnetoresistive Lai-^CaajMnOs system (Qiu et al, 2005; Bozin 
et al, 2006) or ferroelectric lead-based perovskites (Dmowski et al, 2000; Juhas 
et al, 2004). The PDF contains this local information due to the inclusion of dif- 
fuse scattering intensities in the Fourier transform and it is possible to focus the 
modeling on a specific length-scale when searching for matching structure models, 
allowing in principle structure solutions of local, intermediate, and long-range order 
to be obtained separately. 

The procedure assumes a periodic system with known lattice parameters and stoi- 
chiomctry, otherwise there is no information on location or symmetry of the atom sites 
in the unit cell. To solve the unit cell structure the technique constructs a series of trial 
clusters using the PDF-extracted distances. The tested structures are created with a 
direct use of distance information in the experimental data giving it significantly bet- 
ter performance than procedures that search by random structure updates such as 
Monte Carlo based minimization schemes (Juhas et al, 2006; Juhas et al, 2008). 

2. Experimental procedures 

The extended Liga procedure has been tested with experimental x-ray PDFs collected 
from inorganic test materials. Powder samples of Ag, BaTiOa, C-graphite, CaTi03, 
CdSe, Ce02, NaCl, Ni, PbS, PbTe, Si, SrTiOa, Ti02 (rutile), Zn, ZnS (sphalerite) 
and ZnS (wurtzite) were obtained from commercial suppliers. Samples were ground 
in agate mortar to decrease their crystallite size and improve powder averaging. The 
experimental PDFs were measured using synchrotron x-ray diffraction at the 6ID-D 
beamline of the Advanced Photon Source, Argonne National Laboratory using the x- 
ray energies of 87 and 98 keV. The samples were mounted using a thin kapton tape in 
a circular, 10 mm hole of a 1 mm thick flat plate holder, which was positioned in trans- 
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mission geometry with respect to the beam. The x-ray data were measured using the 
"Rapid Acquisition" (RA-PDF) setup, where the diffracted intensities were scanned 
by a MAR345 image plate detector, placed about 20 cm behind the sample (Chupas 
et al., 2003). All measurements were performed at room temperature. The raw detec- 
tor images were integrated using the Fit2D program (Hammersley, 1998) to reduce 
them to a standard intensity vs. 26 powder data. The integrated data were then 
converted by the PDFgetX2 program (Qiu et al., 2004) to experimental PDFs. The 
conversion to PDF was conducted with corrections for Compton scattering, polariza- 
tion and fluorescence effect, as available in the PDFgetX2 program. The maximum 
value of the scattering wavevcctor Qmax rang cd from 19 A"^ to 29 1"^ based on a 
visual inspection of the noise in the F{Q) = Q[S{Q) — 1] curves. 

The PDF function G{r) was obtained by a Fourier transformation of F{Q), 

2 /-Qmax 

G{r) = - F{Q)smQrdQ, (1) 

71" JQmin 

and provided a scaled measure of finding a pair of atoms separated by distance r 

Gir) = Y.i+i filAr - n,) - 47rpor. (2) 

The G{r) function has a convenient property that its peak amplitudes and standard 
deviations remain essentially constant with r and is thus suitable for curve fitting. A 
detailed discussions of the PDF theory, data acquisition and applications for structure 
analysis can be found in (Egami & Billinge, 2003; Farrow & Billinge, 2009). 

3. Methods 

The structure solution procedure was carried out in three separate steps, as described 
in the sections below. The first step consists of peak search and profile fitting in the 
experimental PDF to identify prominent inter-atomic distances up to a cutoff distance 
dcut ■ We have developed an automated peak extraction method which eases this task. 
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In the second step these distances are used as inputs for the Liga algorithm, which 
searches for unit cell positions that give structure with the best match in pair lengths. 
If the sample has several chemical species, a final "coloring" step is necessary to assign 
proper atom species to the unit cell sites. This can be done by making use of PDF 
peak amplitude information. However, we have found that coloring can be also solved 
by optimizing the overlap of the empirical atom radii at the neighboring sites, which 
is simpler to implement and works with greater reliability. 

To verify the quality and uniqueness of the structure, the Liga algorithm has been 
run for each sample multiple (at least 10) times with the same inputs, but different 
seeds of the random number generator. For most samples the resulting structures 
were all equivalent, but sometimes the program gave several geometries with similar 
agreement to the PDF-extracted pair distances. In all these cases the correct structure 
could be resolved in the coloring step, where it displayed significantly lower atom radii 
overlap and converged to known structure solution. A small number of structures 
would not solve by this process and the reasons for failure are discussed below. 

3.1. Extraction of pair distances from the experimental PDF 

In the PDF frequent pair distances generate sharp peaks in the measured G{r) curve 
with amplitudes following Equation (2). The peaks are broadened to approximately 
Gaussian shape that reflects atom thermal vibrations and limited experimental reso- 
lution. Additional broadening and oscillations are introduced to the PDF due to the 
maximum wavevector Qmax that can be achieved in the measurement. This cutoff in 
Qmax in effect convolutes ideal peak profiles with a sine function sin((5max thus 
creating satellite termination ripples. 

Recovering the underlying peaks from the PDF is not trivial. The experimental 
curve can have false peaks due to termination ripples. Nearby peaks can overlap and 
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produce complicated profiles that are difficult to decompose. To simplify the process of 
extracting inter-atomic distances we have developed an automated method for peak 
fitting that adds peak profiles to fit the data to some user-defined tolerance, while 
using as few peaks as possible to avoid over-fitting. This method grows peak-like 
clusters of data points while fitting one or more model peaks to each cluster. Adjacent 
clusters iteratively combine until there is a single cluster with a model that fits the 
entire data set. This allows a steady growth in model complexity by progressively 
refining earlier and less accurate models. Furthermore, most adjustable parameters 
can be estimated, in principle, from experimental knowns. A full description of the 
peak extraction method will be presented in a future paper. 

The present work uses the simplest model for peaks, fitting the G{r) data with 
Gaussian peaks over r and using an assumed value of po- This model ignores the 
effect of termination ripples, but for our data the spurious peaks due to these ripples 
were usually identifiable by their small size. Furthermore, the Liga algorithm is not 
required to use every distance it is given, and should exhibit a limited tolerance of 
faulty distances. The peak fitting procedure returns positions, widths and integrated 
areas of the extracted peaks, of which only the peak positions were used for structure 
determination. 

The peak extraction procedure was implemented in Mathematica 6 and tested on 
both the experimental and simulated data. A typical runtime was about 5 minutes. 
Since the structures are known we can compare the results of the peak extraction 
with the expected results. For both experimental and simulated PDFs of the tested 
structures these compared qualitatively well to the ideal distances up to ~10-15 A, 
including accurate identification of some obscured peaks. Past that range the number 
of distinct, but very close, distances in the actual structure is so great that reliable peak 
extraction is much more difficult. For this reason we only performed peak extraction up 
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to 10 A before running the trials described in Section 4. Apart from removing peaks 
below a noise threshold in order to filter termination ripples out, and one difficult 
peak in the graphite data, all distances used in the structure solution trials below 
come directly from the peak extraction method. 

3.2. Unit cell reconstruction using the Liga algorithm 

In the second step the Liga algorithm searches for the atom positions in the unit cell 
that make the best agreement to the extracted pair distances. The quality of distance 
match is expressed by cost Cd defined as a mean square difference between observed 
and modeled pair distances. 

Cd = ^ J2dk«icut i^k,near - dkf (3) 

The index k goes over all pair distances dfe in the model that are shorter than the 
cutoff length dcut and compares them with the nearest observed distance tk^neari while 
P is the number of model distances. This cost definition considers only distance values 
as extracted from the PDF peak positions, and ignores their relative occurrences. For 
multi-component systems there is in fact no straightforward way of extracting distance 
multiplicities, because it is not known what atom pairs are present in each PDF 
peak. Nevertheless, the cost definition still imposes strict requirements on the model 
structure, as displayed in Fig. 1. A site in the unit cell must be at a good, matching 
distance not only from all other cell sites, but also from all of their translational images 
within the cutoff radius. 
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Fig. 1. Schematic calculation of the distance cost C^. To achieve low a unit cell site 
needs to be at a correct distance from other cell sites and from their translational 
images. 

To find an optimum atom position in the unit cell structure the Liga algorithm 
uses input pair distances in an iterative build-up and disassembly of partial struc- 
tures (Juhas et al, 2006). The procedure maintains a pool of partial unit cell structures 
at each possible size from a single atom up to a complete unit cell. These "candidate 
clusters" are assigned to "divisions" according to the number of sites they contain, 
therefore there are as many divisions as is the number of atoms in a complete unit cell. 
The structures at each division compete against each other in a stochastic process, 
where the probability of winning equals the reciprocal distance cost Cd, and a win is 
thus more likely for low-cost structures. The winning cluster is selected for "promo- 
tion," where it adds one or more atoms to the structure and thus advances to a higher 
division. At the new division a poorly performing high-cost candidate is "relegated" 
to the original division of the promoted structure, thus keeping the total number of 
structures at each division constant. The relegation is accomplished by removing cell 
sites that have the largest contributions to the total cost of the structure. Both pro- 
motion and relegation steps are followed by downhill relaxation of the worst site, i.e., 
the site with the largest share of the total cost C^- The process of promotion and 

lUCr macros version 2.1.1: 2009/05/11 



9 

relegation is performed at every division in a "season" of competitions. These seasons 
are repeated many times until a full sized structure attains sufficiently low cost or 
until a user-specified time limit. A complete description of the Liga algorithm details 
can be found in (Juhas et al, 2008). 

3.3. Atom assignment 

The Liga algorithm used in the structure solution step has no notion of chemical 
species and therefore returns only coordinates of the atom sites in the unit cell. For 
a multi-component system an additional step, dubbed coloring, is necessary to assign 
chemical elements over known cell sites. To assess the quality of different assignments 
we have tested two definitions for a cost of a particular coloring. The first method 
uses a weighted residuum, Ryj, from a least-squares PDF refinement to the input 
PDF data (Egami & Billinge, 2003). The PDF refinement was performed with a fully 
automated PDFfit2 script, where the atom positions were all fixed and only the atomic 
displacement factors, PDF scale factor and Q-resolution damping factor were allowed 
to vary. The second procedure defines coloring cost Cc as an average overlap of the 
empirical atomic radii, so that 

1 2 
Cc=J^ Ed,<r,,i+r,,2 irk,i + rk,2 - dk) (4) 

The index k runs over all atom pairs considering periodic boundary conditions, rk,i 
and rfc^2 are the empirical radii values of the first and second atom in the pair fc, and 
N is the number of atoms in the unit cell. 

Considering an N atom structure with s different atom species, the total number of 
possible assignments is given by the multinomial expression N\/{ni\ • • • ?^sO- ■^^'^ 
a 1:1 binary system the number of possible assignments tends to 2^ with increasing 
N. Such exponential growth in possible configurations makes it quickly impossible to 
compare them all in an exhaustive way. We have therefore employed a simple downhill 
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search, which starts with a random element assignment. The initial coloring cost Cc is 
calculated together with a cost change for every possible swap of two atoms between 
unit cell sites. The site flip that results in the largest decrease of the total coloring cost 
is accepted and all cost differences are evaluated again. The site swap is then repeated 
until a minimum configuration is achieved, where all site flips increase the coloring 
cost. The downhill procedure was verified by repeating it 5 times using different initial 
assignments. In nearly all cases these runs converged to the same atom configurations. 

The downhill procedure was performed using both definitions of the coloring cost. 
For the coloring cost obtained by PDF fitting the procedure was an order of magnitude 
slower and less reliable, as the underlying PDF refinements could converge badly for 
poor atom assignments. The second method, which calculated cost from radii-overlap, 
was considerably faster and more robust. For all tested materials, the overlap-based 
coloring assigned all atoms correctly, when run on correct structure geometry. The 
overlap cost was evaluated using either the covalent radii by (Cordero et al., 2008) 
or the ionic radii from (Shannon, 1976) for more ionic compounds. For some ions the 
Shannon table provides several radii values depending on their coordination number or 
spin state. Although these variants in ionic radii can vary by as much as about 30%, 
the choice of particular radius had no effect on the best assignment for all studied 
structures. 

4. Results 

The experimental x-ray PDFs were acquired from 16 test samples with well known 
crystal structures. To verify that the measured PDF data were consistent with estab- 
lished structure results structure refinements of the known structures were carried 
out using the PDFgui program (Farrow et al., 2007). The PDF fits were done with 
structure data obtained from the Crystallography Open Database (COD) (Grazulis 
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et al., 2009). The structure parameters were all kept constant in the refinements, 
which modified only parameters related to the PDF extraction, such as PDF scale, Q 
resolution dampening envelope and a small rescaling of the lattice parameters. These 
refinements are summarized in Table 1, where low values of the fitting residual 
confirm good agreement between experimental PDFs and expected structure results. 

The PDF datasets were then subjected to the peak search, Liga structure solution 
and coloring procedures as described above. To check the stability of this method, 
several structures were solved using an enlarged periodicity of 1x1x2, 1x2x2 or 
2x2x2 super cells. The lattice parameters used in the Liga crystallography step were 
obtained from the positions of the nearest PDF peaks. In several cases, such as for 
BaTiOa where peak search could not resolve tetragonal splitting, the cell parameters 
were taken from the respective GIF reference, as listed in Table 1. 

The structure solution was considered successful if the found structure displayed 
the same nearest neighbor coordination as its GIF reference and no site was ofi^set by 
more than 0.3 A from its correct position. The solution accuracy was evaluated by 
finding the best overlay of the found structure to the reference GIF data. The optimum 
overlay was obtained by an exhaustive search over all symmetry operations defined in 
the GIF file and over all mappings of solved atom sites to all reference sites containing 
the same element. The overlaid structures were then compared for the differences in 
fractional coordinates and for the root mean square distortion Sr of the solved sites 
from their correct positions. Table 2 shows a summary of these results for all tested 
structures. 
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Table 1. List of measured x-ray PDFs and their fitting residua with respect to established 

structures from the literature. 



sample 


Rw 


GIF reference 


Ag 


0.095 


(Wyckoff, 1963) 


BaTiOa 


0.123 


(Mcgaw, 1962) 


C (graphite) 


0.248 


(Wyckoff, 1963) 


CaTiOs 


0.083 


(Sasaki et at, 1987) 


CdSc 


0.149 


(Wyckoff, 1963) 


CcOa 


0.098 


(Wyckoff, 1963 


NaCl 


0.161 


(Jurgens et al, 2000) 


Ni 


0.109 


(Wyckoff, 1963) 


PbS 


0.085 


(Ramsdell, 1925) 


PbTe 


0.070 


(Wyckoff, 1963) 


Si 


0.085 


(Wyckoff, 1963) 


SrTiOs 


0.143 


(Mitchell et al, 2002) 


Ti02 (rutile) 


0.146 


(Meagher & Lager, 1979) 


Zn 


0.105 


(Wyckoff, 1963) 


ZnS (sphalerite) 


0.102 


(Skinner, 1961) 


ZnS (wurtzite) 


O.I74I 


(Wyckoff, 1963) 



refined as mixture of wurtzite and sphalerite phases 



The procedure converged to a correct structure for 14 out of 16 studied samples and 
failed to find the remaining 2. The convergence was more robust for high-symmetry 
structures, such as Ag (f.c.c), NaCl or ZnS sphalerite, which could be reliably solved 
also in enlarged [222] supercells. For all successful runs the distance cost of the 
Liga-solved structure was comparable to the one from the GIF reference and the 
atom overlap measure Cc was close to zero. ZnS sphalerite shows a notable difi^erence 
between the Cj, values of the solution and its GIF reference, however this was caused 
by using a PDF peak position as a cell parameter for the solved structure. Apparently 
the PDF peak extracted at r a; a was slightly off^set with respect to other peaks, 
nevertheless the Liga algorithm still produced atom sites with correct fractional coor- 
dinates. The mean displacement Sr for ZnS is A, because solved structures and GIF 
references were compared using lattice parameters rescaled to their GIF values. 
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Table 2. Summary of tested structure solutions from x-ray PDF data 



sample atoms 


cost Cd (0.01 A2) 


cost Ce (A2) 


deviation of coordinates 




(supercell) 




Liga 


GIF 


Liga 


GIF 


Sx 


Sy 


Sz 


Sr (A) 


successful solutions 




















Ag [111] 


4 


0.0232 


0.136 





0.001 














Ag [222] 


32 


0.0097 


0.136 





0.001 


0.00025 


0.00024 


0.00003 


0.0014 


BaTiOs [111] 


5 


0.370 


0.394 


0.040 


0.042 


0.0057 


0.0066 


0.014 


0.064 


BaTiOg [112] 


10 


0.392 


0.394 


0.058 


0.042 


0.00023 


0.039 


0.018 


0.16 


C graphite [111] 


4 


0.396 


0.574 


0.010 


0.016 


0.0029 


0.0029 


0.036 


0.14 


C graphite [221] 


16 


0.420 


0.574 


0.010 


0.016 


0.0086 


0.0065 


0.036 


0.15 


CdSc [111] 


4 


0.107 


0.138 





0.001 








0.0055 


0.027 


CdSe [221] 


16 


0.0856 


0.138 





0.001 


0.00010 


0.00013 


0.0057 


0.028 


Ce02 [HI] 


12 


0.515 


0.554 




















NaCl [111] 


8 


1.75 


1.71 




















NaCl [222] 


64 


1.20 


1.71 








0.00031 


0.00031 


0.00035 


0.0032 


Ni [111] 


4 


0.0024 


0.0024 




















Ni [222] 


32 


0.0025 


0.0024 








0.00015 


0.00013 


0.00013 


0.0008 


PbS [111] 


8 


0.0125 


0.0104 


0.010 


0.011 














PbS [222] 


64 


0.0140 


0.0104 


0.010 


0.011 


0.00005 


0.00004 


0.00005 


0.0005 


PbTe [111] 


8 


0.0024 


0.0127 


0.097 


0.090 














PbTe [222] 


64 


0.0022 


0.0127 


0.097 


0.090 


0.00011 


0.00011 


0.00008 


0.0011 


Si [111] 


8 


0.0045 


0.0045 




















Si [222] 


64 


0.0048 


0.0045 








0.00010 


0.00009 


0.00008 


0.0009 


SrTiOs [111] 


5 


0.437 


0.437 


0.002 


0.002 














Zn [111] 


2 


0.495 


0.470 














0.027 


0.095 


Zn [222] 


16 


0.564 


0.470 








0.00010 


0.00006 


0.020 


0.080 


ZnS sphalerite [111] 


8 


0.150 


0.0647 




















ZnS sphalerite [222] 


64 


0.160 


0.0647 








0.00029 


0.00033 


0.00031 


0.0028 


ZnS wurtzite [111] 


4 


0.141 


0.152 














0.0038 


0.017 


ZnS wurtzite [221] 


16 


0.165 


0.152 








0.00003 


0.00002 


0.0039 


0.017 


failed solutions 




















CaTiOs [111] 


20 


0.4967 


0.902 


0.52 


0.072 


0.16 


0.14 


0.17 


1.6 


Ti02 rutile [111] 


6 


0.5358 


0.758 


0.40 


0.009 


0.081 


0.24 


0.00004 


0.94 



Cd, Cc - distance and atom overlap cost as defined in equations (3), (4) 

Sx, Sy, Sz - standard deviation in fractional coordinates normalized to a simple [111] cell 
Sr (A) - root mean square displacement of the solved sites from the reference GIF positions 



The structure determination did not work for 2 lower-symmetry samples of CaTiOa 
and Ti02 rutile. In both of these cases, the simulated structure showed significantly 
lower distance cost Cd while its atom overlap Cc was an order of magnitude higher 
than for the correct structure and clearly indicated an unphysical result. Such results 
were caused by a poor quality of the extracted distances, which contained significant 
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errors and omissions with respect to an ideal distance list. The peak search and dis- 
tance extraction is more difficult for lower symmetry structures, because their pair 
distances are more spread and produce small features that can be below the technique 
resolution. Because of poor distance data, the Liga algorithm converged to incorrect 
geometries that actually displayed a better match with the input distances. Both 
CaTiOa and Ti02 were easily solved when run with ideal distances calculated from 
the GIF structure. 

The results in Table 2 suggest several ways to extend the method and improve its 
success rate. First, the Liga geometry solution and coloring steps can be performed 
together, in other words the structure coloring step needs to be merged to a chemistry 
aware Liga procedure. Since atom overlap cost Cc is meaningful and can be easily 
evaluated for partial structures, the total cost minimized by the Liga algorithm should 
equal a weighted sum of Cc and distance cost Cd- Such a cost definition would steer 
the Liga algorithm away from faulty structures found for CaTiOa and Ti02 rutile, 
because both of them had huge atom overlaps Cc- Another improvement is to perform 
PDF refinement for a full sized structure and update its cost formula so that the PDF 
fit residuum Rw is used instead of distance cost C^- Such modification would prevent 
the cost advantage for wrong structures due to errors and omissions in the extracted 
distances. The assumption is that the distance data are still good enough to let the 
Liga algorithm construct the correct structure in one of its many trials. Finally, the 
cost definition for partial structures can be enhanced with other structural criteria such 
as bond valence sums (BVS) agreement (Brese & Okeeffe, 1991; Norberg et al, 2009). 
Bond valence sums are not well determined for incomplete intermediate structures and 
thus cannot fully match their expected values. However, BVS are always increasing, 
therefore a BVS of some ion that is significantly larger than its expected value is a 
clear sign of such a partial structure's poor quality. 
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5. Conclusions 

We have demonstrated the Liga algorithm for structure determination from PDF can 
be extended from its original scope of single-element non-periodic molecules (Juhas 
et al, 2006; Juhas et al, 2008) to multi-component crystalline systems. The proce- 
dure assumes known lattice parameters and it solves structure geometry by optimizing 
pair distances to match the PDF extracted values, while the chemical assignment is 
obtained from minimization of the atomic radii overlap. The procedure was tested on 
x-ray PDF data from 16 test samples, of which in 14 cases it gave the correct struc- 
ture solution. These arc promising results, considering the technique is at a prototype 
stage and will be further developed to improve its ease of use and rate of conver- 
gence. The procedure can be easily amended by a final PDF refinement step. Such an 
implementation could significantly reduce the overhead in PDF analysis of crystalline 
materials, because its most difficult step, a design of suitable structure model, would 
become fully automated. 
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